Astronomy & Astrophysics manuscript no. 14210 © ESO 2010 

July 19, 2010 



o 

(N 

6 



o 

O 

o 



Radiative transfer with scattering for domain-decomposed 3D MHD 
simulations of cool stellar atmospheres 

Numerical methods and application to the quiet, non-magnetic, surface of a 

solar-type star 

W. Hayek>'23, M. Asplund\ M. Carlsson^ R. Trampedach'^'^, R. Collet\ B. V. Gudiksen^ V. H. Hansteen^ and J. 

Leenaarts^ 

' Max Planck Institut ftir Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany 

e-mail: [hayek , asplund , remo] @mpa-garching . mpg . de 
- Research School of Astronomy & Astrophysics, Cotter Road, Weston Creek 26 11, Australia 

^ Institute of Theoretical Astrophysics, University of Oslo, RO. Box 1029 Blindern, N-0315 Oslo, Norway 

e-mail: [mats . carlsson , b . v . gudiksen , viggo . hansteen] ©astro . uio. no 
* JILA, University of Colorado, 440 UCB, Boulder, CO 80309-0440, U.S.A. 

e-mail: trampedaOlcd . Colorado . edu 
^ Astronomical Institute, Utrecht University, Postbus 80 000, NL3508 TA Utrecht, The Netherlands 

e-mail: j.leenaarts@uu.nl 

My 19, 2010 

ABSTRACT 

Aims. We present the implementation of a radiative transfer solver with coherent scattering in the new BIFROST code for radiative 
magneto-hydrodynamical (MHD) simulations of stellar surface convection. The code is fully parallelized using MPl domain decom- 
position, which allows for large grid sizes and improved resolution of hydrodynamical structures. We apply the code to simulate the 
surface granulation in a solar-type star, ignoring magnetic fields, and investigate the importance of coherent scattering for the atmo- 
spheric structure. 

Methods. A scattering term is added to the radiative transfer equation, requiring an iterative computation of the radiation field. We 
use a short-characteristics-based Gauss-Seidel acceleration scheme to compute radiative flux divergences for the energy equation. The 
effects of coherent scattering are tested by comparing the temperature stratification of three 3D time-dependent hydrodynamical at- 
mosphere models of a solar-type star: without scattering, with continuum scattering only, and with both continuum and line scattering. 
Results. We show that continuum scattering does not have a significant impact on the photospheric temperature structure for a star like 
the Sun. Including scattering in line-blanketing, however, leads to a decrease of temperatures by about 350 K below log,Q T5000 ^ -4. 
The effect is opposite to that of ID hydrostatic models in radiative equilibrium, where scattering reduces the cooling effect of strong 
LTE lines in the higher layers of the photosphere. Coherent line scattering also changes the temperature distribution in the high 
atmosphere, where we observe stronger fluctuations compared to a treatment of lines as true absorbers. 

Key words. Radiative transfer - Stars: atmospheres - Sun: atmosphere 



^ . 1. Introduction Assuming a plane-parallel or spherical-symmetric stratification, 

; '~j they include only a rudimentary treatment of convective en- 

;X; The atmospheres of late-type stars form the transition from the ^^gy transport in cool stellar atmosp heres. Subs equent updates 

H opaque convective envelope to the interstellar medium. Hot ris- ^f ^j^^^^ ^^^^^^ (^ g^ j^^^^.^^^ ^996; G ustafsson et al. 2008) ben- 

. >y ing plasma transports heat to the surface, becomes transparent g^j f^^m the largely increased computational power, refining the 

and looses its entropy through radiative cooling. Gravity accel- treatment of the strongly wavel ength-dependent line opacities. 

erates the cooled gas back into the star, carrying kinetic energy ^ewer codes, such as PHOENIX (iHauschildt et al.l[T999l) can also 

inward and forcing the convective flow. By taking over heat include departures from local thermodynamic equilibrium (LTE) 

transport and removing entropy, the radiatiori field therefore in- j^ ^^e radiative transfer computation and the absorber popula- 

directly drives convection d Stem &Nordlund| | 1998|) , making ra- ^^^^^ j^^ ^I^ ^^^^^^ ^^^^ ^^^ ^^jy provided growing insight 

diative and hydrodynamical processes equally important at the ^^^^ j^e physical environment at the surface of cool stars, but 

surface. Magnetic fields have strong impact on the higher atmo- ^^^^ ^j^^ ^^^^^^ ^ ^t^jj^^^^ t^^l ^^ chemical abundance analy- 

sphere and cause local phenomena in the surface granulation, ^^^ ^he wide variety of applications includes studies of galactic 

such as spots and pores. chemical evolution and of the origin of the elements. 

The classical numerical models of cool stellar atmospheres 
in ID focused on a detailed description of radiative trans- The advent of fully dynamic 3D surface convection simu- 

fer, with two promine nt examples being th e MARCS cod e lations has enabled a much more realistic treatment of the hy- 

dGustafsson et al.l Il975h and the ATLAS code (lKuruczlll979h . drodynamical plasma flow, deepening our understanding of con- 
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vection and eliminating the need for microturbulent and macro- 
turbulent broadening in line formation computations (see, e.g., 
[Nordlund et al. 2009). The 3D models are capable of accurately 
reproducing the surface structure of the observed solar gran- 
ulation vyith their strong ly inhomogeneous surface intensities 
JStein & Nordlundl[T998b . The velocity fields predicted by the 
3D simulations lead to a close match with both the observed 
spectral line bisectors and the broadening of th eir profiles in the 
atmospheres of different stars (e.g.,iP ravins & Nordlundl 1990t 
Asplund et aP l2000t lAllende Prietoet aL i2002t iRamirez et alJ 



2009). Recently, impressive agreement between a new synthetic 
3D model and solar observations has been found in a detailed 
comparison of spectral line shifts, equivalent widt hs and center- 
to-limb variations for normalized line profiles dPereira et alj 
l2009alfa) . In essentially all cases, this 3D model reproduced the 
observations with an accuracy that is compara ble to the semi- 
empirical model of iHolweger & Muelleij (1 19741) . which is tradi- 
tionally used in spectroscopy of the solar photosphere. 

The accuracy of the treatment of radiation in 3D, however, 
is still strongly limited by the available computational power 
Radiative transfer easily becomes the most computationally ex- 
pensive part of a simulation, since the equations must be solved 
for a considerably larger set of transport directions compared to 
hydrodynamics, and non-grey opacities must be accounted for 
in realistic simulations. Most of the cuiTently existing 3D ra- 
diative (M)HD codes therefore assume LTE and capture the at- 
mospheric height dependence of conti nuum and line opacities 
using the opacity binning method (e. g.. iNordlundl 1 982t iLudwig 
Il992h : the problem of computing the monochromatic radiation 
field for a larger number of wavelengths is reduced to the nu- 
merical solutio n of the ra diative transfer equation for typically 
5 opacity bins. ISkartlieiJ (|2000) extended the opacity binning 
method to include coherent scattering, and showed its impor- 
tance in the solar chromosphere using a 3D radiative transfer 
solver for parallel shared-memory architectures. 

Modern large-scale computer clusters use distributed mem- 
ory architectures to handle the growing complexity of scien- 
tific simulations, allowing, e.g., self-consistent MHD models of 
the solar chromosphere, tra nsition region and corona (Hansteen 
l2004t iHanstee n et al. 2007) or detailed hydrodynamical models 
of giant stars CCollet et al...2007.) . We present a new fully MPI- 
parallelized radiative transfer solver with coherent scattering for 
the new BIFROST code for time-dependent 3D MHD simulations 
of cool stellar atmospheres (Gudiksen et al., in preparation). 

In Sect.|2]and[3] we discuss the physics of the radiative trans- 
fer model and its implementation in the MHD code. Section |4] 
describes the most important continuous and line opacity sources 
that we include in our simulations. Section|5]describes the appli- 
cation of the BIFROST code to model the atmosphere of a solar- 
type star using radiative transfer calculations with scattering, and 
discusses the effects on the temperature structure. 



2. Radiative transfer with scattering and the 
radiative flux divergence 

2.1. The radiative transfer equation 

Hydrodynamical simulations of cool stellar atmospheres need to 
cover several pressure scale heights above and below the optical 
surface to minimize the effect of the boundaries on the granu- 
lation flow. The exponential density stratification causes the op- 
tical depth of the plasma to span about 15 orders of magnitude 
from the highest to the lowest layers of the simulation. The ra- 
diative transfer problem must therefore be solved in very differ- 



ent physical environments: in the extremely optically thick dif- 
fusion region at the bottom of the simulation box, all photons 
are thermalized. At the top, the atmosphere is mostly optically 
thin and mainly photons in the strongest lines interact with the 
gas. For the bulk of the photons, the transition between these 
two domains is rapid; it is confined to a thin layer which appears 
coiTugated due to the different geometrical depth variation of 
opacities in upflows and downflows (Stein & Nordlund 1998). 

Radiative transfer is, in general, a time-dependent process, 
which needs to be treated simultaneously with the hydrodynam- 
ics. However, the timescale of photon propagation over a mean 
free path length, t^ - {cxa)\ where ^^ is the monochromatic 
opacity and c is the speed of light, is orders of magnitude shorter 
than any hydrodynamical timescale. Radiative transfer therefore 
decouples from the hydrodynamics and is well approximated by 
a time-independent problem, described by a radiative transfer 
equation for the monochromatic specific intensity Ia(x, h) in di- 
rection h: 



h • VIa(x, n) = -xa(x)Ia(x, n) + Ja(x, n). 



(1) 



where /'i denotes the local emission at wavelength A (see, e.g., 
Mihalas 1978; Rutt'enl l2003h . The left-hand side of Eq. dB is de- 
fined in the rest frame of the model atmosphere. The source and 
sink terms of the right-hand side are naturally described in the 
co-moving frame of the flowing gas. The consequent Doppler 
shifts are difficult to treat in 3D time-dependent simulations due 
to restrictions in computational power, requiring us to compute 
wavelength-integrated quantities in the opacity binning approx- 
imation (see below). We therefore assume a static medium, ne- 
glecting possible influences of the velocity field. 

The extinction of photons is described, as customary, 
through the absorption coefficient k,i and the scattering coeffi- 
cient o-,i, which combine to the gas opacity. 



Xa(x) = ka(x) + o-a(x). 



(2) 



and give rise to the definition of the photon destruction probabil- 
ity 



q(-k) 



ka(x) 



ka{x) + o-a{x) 



(3) 



Recasting the optical path ds - h • dx along a ray in direction 
h into the optical depth cIta - XAds along that direction, gives 
Eq. ([TJ the form 



dh 
dr. 



■(ta) = -Ia(ta) + S a{ta), 



(4) 



with the source function S a = JaIxa- For the numerical compu- 
tation, we employ the formal solution 



hiTA) = /i(Tu,i)e 



(T,l-T„„,) ^ I 



S A(t)e'-'''dt, 



(5) 



where Ia{tu,a) is the incident intensity at the upwind end of the 
ray at optical depth t^j < ta- 

The source function S a at optical depth ta in direction n in- 
cludes local thermal radiation from the gas and coherent scatter- 
ing of photons: 



Sa^ 



Q-A 
47r;r. 



■ I 0(«,«')/.i,n't/Q' + — B,i = (1 -e,07^ + e^B,i, 
Js^ Xi 



(6) 



where scattered radiation from direction h' contributes with 
weight (p in the integral over the unit sphere S ^, Ba denotes the 
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Planck function, and J^ is the mean intensity. The second equal- 
ity holds for isotropic angular redistribution of radiation {(p = 1). 
For 6) < 1, the source function depends on J,i and thus, through 
the non-locality of radiative transfer, on radiation processes in 
the entire simulation domain. This turns Eq. (|4]l from an ordi- 
nary differential equation into an integro-differential equation. 

CuiTent limitations of available computing resources require 
the assumption of isotropic coherent scattering. Continuum 
processes in cool stellar atmospheres and very strong lines 
fulfill this restriction in very good or reasonable approxima- 
tion, respectively, due to their weak wavelength dependence. 
Intermediate and weak lines are more accurately treated in com- 
plete spectral redistribution. 

2.2. The radiative flux divergence and the wavelength 
integral 

Absorption and thermal emission of radiation couples the stel- 
lar plasma with the radiation field through the transfer of heat. 
Photon energies in cool stars are too small to exert a significant 
force on the fluid compared to the gas pressure and gravity; the 
coupling is therefore sufficiently described by adding a radiative 
heating term 2,ad to the energy equation. 

Evaluating the first moment of Eq. ([T]i and using the above 
definitions yields 



- V . F, = AnxA (J A -Sa)^ 4neAXA (J a - Ba) , 



(7) 



where V • Fa is the local monochromatic radiative flux diver- 
gence. The second equality holds in the case of the coherent 
scattering source function (Eq. (|6]l). The scattering term does 
not contribute to heat exchange by definition, reducing radiative 
heating and cooling by a factor of q compared to the case where 
Sa = Ba- 

Integrating the monochromatic flux divergence in Eq. O 
over the whole wavelength spectrum of the star yields the local 
heating rate grad: 



G, 



rad 



V • F^dA. 



(8) 



In the optically thick regime, where radiative transfer is diffu- 
sive, this integral may be simplified with good accuracy by as- 
suming the Rosseland mean opacity in the so-called gray ap- 
proximation. However, gray opacities are not sufficient for a re- 
alistic treatment of the height-dependent line-blanketing above 
the surface, where the atmospheric structure is very sensitive 
to the radiation field. Atomic and molecular lines are important 
opacity sources in this region, changing the radiative heating and 
cooling c ompared to the sim plified case of a gray atmosphere 
(see, e.g.. lVogler e t al."2004l for a detailed discussion). The cur- 
rent version of the MARCS I D atinosph ere code uses the opacity 
sampling technique ( Pevtre niannl 1 974 ). which approximates the 
spectrum through statistical sampling at ~ 100000 wavelength 
points. This resolution is sufficient to capture continuum absorp- 
tion and line-blanketing without bias, at least in the lower parts 
of the atmosphere where the spectral distribution of absorbers is 
sufficiently widespread. Stellar atmosphere models in 3D do not 
allow for such a detailed treatment yet, since a single formal so- 
lution is many orders of magnitude more expensive to compute: 
the radiative transfer equation in our radiation-hydrodynamical 
model (Sect. 15.11 ) is solved for 240 x 240 columns and 24 trans- 
port directions, which is equivalent to ~ 10^ ID calculations for 
each time step. 



iNordlundl (ll982h . lLudwi"gl(ll992h and lSkartlienl(l2000h have 
described opacity binning techniques, where wavelength inte- 
gration is performed over subsets of the spectral range before 
the solution of the radiative transfer equation, and the radiation 
field is computed for only a few mean opacities instead of the full 
spectrum. We will give a brief d escription of the technique in the 
following; see lSkartlien I (2000) for a more detailed discussion. 

Integrating the radiative transfer equation (Eq. ([TJ) over 
wavelength leads to the definition of a mean opacity, mean scat- 
tering coefficient and mean absorption coefficient; 



X 



or' 



JxA^AdA 
JiAdA 

jo- A J Ad A 
jJAdA 

j KABAdA 

jBAdA ' 



(9) 
(10) 
(11) 



The intensity-weighted mean opacity ;(f^ and the mean-intensity- 
weighted mean scattering coefficient o-' depend on the unknown 
radiation field Ia and its angular average Ja, which must be es- 
timated: we use ID radiative transfer calculations on the mean 
stratification of the atmosphere (see Sect. 15.2b . which yield ap- 
proximations for;^f' X x'-^^ and o-' x; o-'-^^. 

These three mean coefficients represent absorption, scatter- 
ing and thermal emission of photons with good accuracy where 
the stellar atmosphere is optically thin across the spectrum. 
However, x''^^ does not ensure a correct total radiative energy 
flux at optical depths r » 1 where radiative transfer is diffusive. 
It needs to be replaced by the Rosseland mean opacity, defined 
as the weighted harmonic mean 



X 



J(dBA/ds)dA 
J{l/XA)idBA/ds)dA' 



(12) 



We consequently use a r-weighted sum of the two quantities 
^j.iD ajjjj^R xjje geometrical depth of the transition between the 
two regimes near r ^ 1 varies quickly with wavelength where 
spectral lines are present, and it is not sufficient to consider 
only a single pair of mean opacities x^'^^ and x^- The opacity 
binning method therefore defines several opacity groups, where 
each member reaches unit optical depth (ta - 1) at a similar 
geometrical depth. The integrals in Eq. (|9j - Eq. (fT2l i are then 
evaluated only for a set of member wavelengths {/I,} in each bin 
/, which does not have to be continuous. 

Depending on the height range of the stellar atmosphere 
model and the wavelength selection method, it turns out that 
about 5 such opacity bins are enough to capture the essence of 
the line-blanketing and continuum opacity and to obtain a real- 
istic temperature structure ( Vogler et al. 2004). M ore recent at- 
mosphere models have been extended to 12 bins (ICaffau et al.l 
2008.) . For the simulations presented in this work, we compute 
radiative transfer with 12 bins, where wavelengths are sorted not 
only by the geometrical height of the monochromatic optical sur- 
face, but also by wavelength, separating opacities in the UV, vi- 
sual and infrared bands (Trampedach et al., in preparation). 

It is difficult to assess the quality of the opacity binning 
method in realistic 3D simulations: deviations of the resulting 
radiative heating rates grad from an accurate monochromatic so- 
lution have a height-dependent impact on the temperature struc- 
ture (see Sect. (|5]l), making the long-term behavior of the simu- 
lation hard to predict. The agreement of 3D model atmospheres 
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with various observational tests indicates that opacity binning 
still yields a reasonable estimate for the line-blanketing. 

3. The numerical implementation 

The large variety of radiative transfer models for astrophysical 
problems inspired the development of very different analytical 
and numerical methods to obtain the radiation field (see, e.g., 
IWehrse & KalkofenI (l2006h for an overview). For our given prob- 
lem of computing radiative heating rates as the flux divergence 
-V • F of a time-independent radiation field in 3D, the direct so- 
lution of Eq. (HJl yields accurate results with efficient numerical 
schemes. 

Characteristics methods, which solve the transfer problem 
along a discrete set of light rays to capture the anisotropy of 
the radiation field in the optically thin at mosphere, are a pop- 
ular choice in s tellar atmosphere models. iNordlundl (Il982h and 
ISkartliedj 20001) u se Feautrier-type differential radiative transfer 
solvers (lFeautrieii ll964) for solving Eq. (|4| on long characteris- 
tics. They span across the entire simulation domain, which is an 
obstacle for a domain-decom posed parall e lizatio n of the MHD 
code (see Sect. [HI below ). iBrulset al.l (Il999l) . IVogler et alj 
(l2005h and LMuthsam et al.' ('2010') employ the short charac- 
teristics meth od (Mihalas et al. 1978; Olson & Kunasz 1987; 
iKunasz & Aueii ll988). where the radiative transfer equation is 
solved on characteristics which only extend to the adjacent up- 
wind and downwind grid layers. This method is required by our 
choice of iteration technique for an efficient solution of the scat- 
tering problem. 



3.1. Short characteristics 

The short characteristics method employs the formal solu- 
tion (Eq. Q) of the monochromatic radiative transfer equation 
(Eq. (|4]i) to compute the radiation field at the center of a three- 
point ray for a known source function S a- The discretization is 
performed by interpolating the source function for a given wave- 
length A (or bin number) along the ra y using a s econd-order 



Bezier curve (see, e.g., the discussion in lAueill2003h 

S(t) = (1 - tfSu + t^So + 2t(l - t)Sc 



(13) 



where Su and So are the upwind and local source functions and 
t - {t - Tu)/(to - Tu) is the curve parameter. Sc is a control 
point, which shapes the interpolating curve by restricting it to 
the convex hull laid out by Su, Sc and Sq. This characteristic 
of Bezier curves may be exploited to detect and suppress over- 
shoots, which destabilize the numerical solution at places in the 
atmosphere where strong opacity and temperature gradients oc- 
cur. Inserting Eq. (fT3l l into the formal solution (Eq. (|5]l), evalu- 
ating the integral and reordering the terms yields the discretized 
expression 



I(t) = /(Tje 



-(T-T„) 



+ >Pu5u+^o5o+*I'd5d. 



(14) 



The shape of the three interpolation coefficients ^'^i, ^'() and *Pd 
for the upwind, center and downwind source functions depends 
on the control point Sc- Choosing 



TO' 



^S' 



(15) 



where 5^ is the centered derivative on the three-point stencil 
(Su,S(),Sd), yields second-order interpolation. It is used where 
no overshoots happen and correctly reproduces the diffusion ap- 
proximation at optical depths t > 30 (see Appendix |A] for the 



detailed shape of the T coefficients). In the optically thin atmo- 
sphere where r < 10"^, a second-order expansion of the (1 -e"^) 
terms in the "P constants stabilizes the solver, which may there- 
fore be implemented with single precision floating point numer- 
ics throughout the simulation domain. Optical depths At along 
the characteristics are similarly computed using the Bezier inter- 
polation technique. 

The mean intensities J and the components of the flux vector 
F are computed by approximating the zeroth and first moment 
integrals by a quadrature sum over selected ray angles ("method 
of discrete ordinates"). 



■^ ~ 4^ X ^'^*^"''^' ^J ~ X ^i^("i'>("i • "7), 



(16) 



where w,- is the weight of direction «,-. The best choice of quadra- 
ture depends on the expected anisotropy of the radiation field 
and on the quantity that needs to be computed. In our case, the 
components of the flux vector F need to be calculated explic- 
itly, requiring the quadrature be invariant to rot ation by 7r/2 to 
avoid directional bias. Carlson's A4 quadrature (ICarlsonlll963h 
with 3 rays per octant is an appropriate choice and represents the 
anisotropy with good accuracy. 

Short characteristics require knowledge of the upwind inten- 
sities /(tu) for each ray direction h, on which the sweep direction 
for a formal solution therefore depends. Interpolation yields all 
such quantities (Sect. 13.31 ). Shallow rays, that fail to hit the up- 
wind layer within the grid cells, need to be extended and may 
cross several cells, possibly across subdomain boundaries. For 
the first formal solution of a simulation run, a Feautrier-type long 
characteristics solver delivers boundary intensity estimates; in- 
tensities from the previous iteration in the neighbor subdomains 
are used for all subsequent computations. Once /(t) is known 
along two edges of the current layer, the remaining unknown 
intensities may be computed away from the boundary through 
vertical interpolation between the upwind layer and the cuiTent 
layer It is worth noting that some long characteristics codes turn 
transport directions around the vertical axis with every time step 
to avoid numerical artefacts stemming from a fixed set of dis- 
crete ordinates. Such an effect is not observed in our short char- 
acteristics implementation. Moreover, the anisotropy of the ra- 
diation field slows down convergence of an iterative solution in 
optically thin parts when transport directions are turned between 
time steps, since the stored boundary intensities come from the 
previous solution (see Sect. 13.21 for further details). All ray di- 
rections are therefore kept fixed. 

The discretized formal solution (Eq. [T4b in the simulation 
domain and averaging of the radiation field over solid angle will 
be abbreviated in the following using the A operator, which is 
commonly defined through 



J^AS. 



(17) 



A is a linear matrix operator on the source functions S which 
represents the numerical algorithm used to compute the radiation 
field in the code. 



3.2. The Gauss-Seidel scheme and MPI paralielization 



As noted in Sect. 12. II the coherent scattering term turns the trans- 
fer equation into an integro-differential equation for the specific 
intensity /. Using the A operator defined in Eq. (fTTI i. the problem 
may be rewritten into the matrix equation 



[l-(l-e)A]S =eB, 



(18) 
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Fig. 1. Horizontal mean photon destruction probability e for three bins representing continuum and weak lines (1), intermediate 
lines (2) and strong lines (3) in the UV, plotted as a function of geometrical depth (left) and optical depth in the respective bin (right, 
the average is taken over surfaces with the same vertical optical depth). The dotted line marks the stellar surface (left) and unit 
optical depth in each bin (right). 



with the identity matrix 1. The expression represents a very large 
system of linear equations. Its direct solution in 3D through in- 
version of the operator on the left-hand side is far too com- 
plex and numerically unstable in some cases, to be of practi- 
cal use. Most solvers therefore apply an iteration scheme, the 
choice of which depends on the structure of the A operator 
matrix. Approximate Lambda Iterati on (ALI, soni etimes also 
called Accelerated Lambda Iteration, ICannon|[l973h is a popu- 
lar method to obtain a good approximation of the radiation field 
with fast convergence. J is computed through a formal solution 
and used to correct the source function. Rather than just insert- 
ing 7 in 5 , which leads to very slow convergence (or no con- 
vergence at all), a largely simplified approximate operator A* is 
used to compute correction values A5 at low cost, speeding up 
convergence tremendously. 

We employ die Gauss- Seidel scheme 

(iTruiillo Bueno & Fabiani Bendichd Il995h . an ALI method 
that combines the formal solution and correction steps. It 
mimics a tridiagonal A* operator, but the scheme does not re- 
quire the expensive construction of the matrix. Source function 
corrections at the grid point / are obtained during a solver sweep 
from the expression 



lOld/new 



A5,= 



l-(l-e,)A,7 



old 



(19) 



yo /new jg jj^g radiation field that includes the coiTections in the 
upwind part of the simulation domain, which have already been 
computed during the current sweep. The dependence of A5 , on 
/, in each layer for immediate correction of S , during the sweep 
requires employing the short characteristics method. The de- 
nominator contains the diagonal element A„ of the A operator, 
which may be computed using Eq. (fT4b and therefore reduces 
to a sum of "P constants. Source function coiTections may be 
applied during both upsweeps and downsweeps for faster con- 
vergence. 

We tested our radiative transfer code by compar- 
ing the numerical results with an analytical solution for 
the case of an isothermal ID atmosphere with constant 
photon destruction probability e (see the discussion in 
iTruiillo Bueno & Fabiani BendichdlI995h and found very good 
agreement. 



The radiation solver is parallelized using spatial domain de- 
composition and communication with the MPI library, adopting 
the virtual topology given by the MHD solver of the BIFROST 
code. The grid is decomposed into cuboid subdomains, allowing 
an arbitrary number of divisions on all three spatial axes. While 
this parallelization lends itself to a mixed initial and bound- 
ary value problem found in computational hydrodynamics, it is 
harder to apply in an efficient way to the pure boundary value 
problem of time-independent radiative transfer Concurrent com- 
putation of spectral subdomains (or opacity bins) would provide 
a higher degree of parallelism considering the non-local depen- 
dencies in a monochromatic formal solution of our given coher- 
ent scattering problem, but such an approach would cause severe 
load balancing issues and suffer from node memory limitations 
when applying the code to very large simulations. Spatial do- 
main decomposition may still be combined with spectral domain 
decomposition if radiative transfer needs to be solved for a large 
number of wavelengths. 



iHeinemann et al.l ( l2006l) have presented a domain- 
decomposed method based on a variant of the formal solution 
(Eq. (|5]l) on long characteristics. The solver bypasses the 
problem of missing incident intensities at subdomain bound- 
aries by splitting the local and boundary contributions. While 
their approach efficiently solves the radiative transfer equation 
without scattering, the long characteristics solver would have to 
be combined with a different ALI scheme than Gauss-Seidel. 
An approximate A* operator needs a certain bandwidth around 
its matrix diagonal to achieve good convergence (see, e.g., the 
discussion in Hauschildt & Baron 2006) . It is therefore more 
expensive to construct and invert than the diagonal operator 
used for the Gauss-Seidel scheme. 

Our code iterates the solution, starting with the source func- 
tion and subdomain boundary intensities from the previous hy- 
drodynamical time step, until the maximum relative source func- 
tion coiTection in the domain after the nth iteration is smaller 
than a preset threshold C: 
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When scattering is not included, the maximum relative change 
of mean intensities at the boundary is used instead to test the 
convergence of the radiation field; 
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If too few iterations are performed, the subdomain boundaries 
produce artifacts in the upper parts of the atmosphere, where 
photon mean free paths are comparable to or larger than the sub- 
domain size. In practice, it turns out that a threshold of C ~ 10"'' 
yields good results in either case. 

The convergence speed of an iterative method depends on 
the spectral radius p of the operator with which corrections 
are computed, as the error of the solution after n iterations de- 
creases with p". The spectral radius approaches p ^ \ - e 
for optically thick scattering media (see, e.g., the discussion in 
iTruiillo Bueno & Fabiani Bendichdll995h . Strong scattering at 
high optical depths therefore leads to very poor convergence 
rates of the Gauss-Seidel solver, requiring hundreds of iterations 
in extreme situations. However, this difficulty is mostly allevi- 
ated by using the source function solution from the previous time 
step and the slow evolution of the plasma flow between consec- 
utive time steps, so that the code ideally needs to fully converge 
the solution only once at the beginning. Domain decomposition 
additionally slows down convergence if the photon mean free 
paths cross subdomain boundaries, which is the case at contin- 
uum wavelengths in the thin atmosphere, since the subdomain 
boundary intensities are not initially known. Storing intensities 
from the previous time step again largely circumvents this prob- 
lem, and the actual number of iterations per time step that is 
required during a simulation run depends on how fast the atmo- 
sphere evolves. 

We therefore test the convergence of the solution for arbi- 
trary time steps of our solar-type simulation using 12 opacity 
bins with continuum and line scatterin g (see Sect.|5]), following a 
similar discussion in lSkartlienI (l2000l) . The tests were run at half 
resolution on all axes to facilitate computation on a single core, 
which yields slightly faster convergence. Since the true solution 
S of our radiative transfer problem is unknown, we compare the 
approximate solution after n iterations, S", with an approximate 
solution S°° which we obtained after additional iterations with 
a lower convergence threshold of C ~ 10 '*, assuming S'^ ~ S 
with good accuracy. 

We use three representative opacity bins, which cover weak, 
intermediate and strong opacities in the UV, with different depth- 
dependence of the scattering strengths. The remaining nine bins 
at longer wavelengths behave in a similar way. Figure [T] shows 
horizontal averages of the photon destruction probabilities e for 
each bin in an arbitrary snapshot of our photospheric simulation: 
averages over layers with the same geometrical depth are plotted 
in the left panel, averages over surfaces with the same vertical 
optical depth are plotted in the right panel. 

Figure |2] compares the convergence speed for the radiative 
transfer solution of the sample bins with and without domain 
decomposition, and with different time step lengths. Thick lines 
represent convergence relative to the true solution S °° for each 
bin, thin lines show the convergence relative to the solution ob- 
tained in the previous iteration, which we use as the convergence 
criterion. In normal operation, the solver would stop as soon as 
the thin line of the currently computed opacity bin has crossed 
the dotted horizontal line. We caution that the number of itera- 
tions needed for a solution also depends mildly on the time step- 
ping algorithm, since the choice of method affects the deviation 



of stored boundary intensities and source functions between sub- 
steps of the time integration. We therefore only analyze the be- 
havior for the first extrapolation step of a 3rd order Runge-Kutta 
time stepper. 

The poorer convergence speed caused by scattering at high 
optical depths in bin 3 is evident in all plots (thick dot-dashed 
line), compared to the situation in bin 1, where the photon de- 
struction probability is larger. The small optical path lengths of 
bin 3 reduce the impact of domain decomposition, since the radi- 
ation field is essentially local in most parts of the simulation box. 
Contrary to that, bin 1 suffers most strongly from slower conver- 
gence with increasing number of subdomain divisions, as well as 
from some flip-flopping of AS' . The latter is caused by high-order 
interpolation (see Sect. l3.3l ) and disappears when the solver is set 
to linear interpolation. High order interpolation of upwind inten- 
sities widens the domain of dependence of the short character- 
istics, and the effect is amplified where large path lengths in the 
optically thin regime cross subdomain boundaries. 

Domain decomposition mildly slows down convergence, and 
the accuracy of the solution in bin 3 slightly deteriorates for a 
larger number of subdomains. Longer time steps have the same 
effect on that bin, causing slower convergence towards S °° than 
indicated by the relative coiTections with respect to 5""' (thick 
and thin dot-dashed lines in Fig. |2l). The method devised by 
iSkartlienI ( |2000|) exhibits similar behavior for bins with strong 
scattering lines. 

The effect of such inaccuracies in the numerical solution of 
S and J on the energy transfer between the radiation field and 
the gas are nevertheless small or even vanish in some regions: 
radiative heating is reduced in the atmosphere where coherent 
scattering is important (see Eq. (|7]) and Fig.[T]). Coherent scatter- 
ing also effectively damps the impact of any remaining disconti- 
nuities in the radiation field across subdomain boundaries on the 
flux divergence in the optically thin atmosphere, so that no vis- 
ible artifacts from the domain decomposition remain in the gas 
temperatures. 

Co mpared to the solver proposed by iHeinemann et al.l 
(l2006h . it is clear that our method is not optimal for the case 
without scattering, since several computationally expensive for- 
mal solution and communication steps are required to obtain a 
radiation field that is consistent in the whole domain. It offers 
good performance when scattering is included, which is not con- 
sidered in their method. 

3.3. Interpolation and grid refinement 

At every time step, the hydrodynamical solver updates mass den- 
sities and internal energy densities. These quantities are used to 
look up tabulated opacities, bin-integrated Planck functions and 
photon destruction probabilities at every grid point. In general, 
the characteristics grid needed to represent the anisotropy of the 
radiation field does not coincide with the hydrodynamical mesh, 
requiring the interpolation of ;t', S and the upwind intensities I^ 
during the formal solution. 

The accuracy of this interpolation strongly influences the 
overall accuracy of the solver, and there is a large choic e of pos- 
sible methods (see, e.g., the discussion in lAueiir2003h . Linear 
interpolation is fast and avoids instabilities produced by interpo- 
lation overshoots, but yields poor estimates where the radiation 
field is not well-resolved, e.g. between granules and intergran- 
ular lanes at the optical surface. It also amplifies the numerical 
diffusion effect of short characteristics, where lateral diffusion 
artificially transports radiation away from the beam. 
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Fig. 2. Convergence of the source function for bins 1-3 (see Fig. [TJ during a simulation run without domain decomposition (left 
column), with 2x2x2 decomposition (center column) and 3x3x3 decomposition (right column), and with a time step of 
Af - 0.03 s (upper row). At - 0.04 s (center row) and At = 0.08 s (lower row). Line styles represent the same bins as in Fig. [1] Thin 
lines: relative source function correction AS after n iterations with respect to 5" ' from the previous iteration « - 1. Thick lines: 
relative source function correction AS with respect to the "true" solution 5°°. Dotted lines mark the threshold C beneath which 
convergence is assumed (see text). 



To illustrate this be havior, we repeat the searchlight test of 
iKunasz & Auerl (Il988h . where a rectangular light beam is cast 
through an empty 3D box with a 100^ mesh and zero opac- 
ity. Any diffusion of radiation away from the beam results in 
a broadening of the beam profile at the surface and can only 
stem from the interpolation of unattenuated upwind intensities. 
The light source illuminates the bottom of the 3D box, where it 
initially covers an area of 30^ mesh points; it is slanted with an 
angle of = 28.1° off the vertical and an azimuth of (/> - 45 .0° . 
The upper left panel in Fig. [3] shows the beam profile at the top 
of the 3D box expected from an exact solution of the unattenu- 
ated transfer problem through vacuum; note that the finite reso- 
lution of the surface in the plot leads to a slightly widened pro- 
file. The upper right panel shows the broadening of the beam 
profile caused by 100 consecutive linear interpolations applied 
for the numerical transfer through the box. Although the area- 
integrated intensity is conserved with good accuracy, limited 
by the machine precision, the beam is visibly widened through 
numerical diffusion. The lower left panel in Fig. [3] shows the 
result when using local cubic interpolation for the transport 
problem. The broadening is reduced, but the overshooting cu- 
bic polynomials produce ringing and negative intensities. We 
th erefore use the local cub ic monotonic interpolation scheme 
of iFritsch & Butland (Il984l) . which effectively suppresses over- 



shoots by using weighted harmonic mean derivatives, in con- 
secutive ID- ID interpolation on horizontal planes, and local 
quadratic interpolation on vertical cell walls (see AppendixiBjfor 
further details). The lower right panel in Fig. |3] shows the result 
from the searchlight test, where the beam profile is conserved to 
a satisfactory degree. Numerical diffusion is reduced and reaches 
a level which renders the computed flux divergences comparable 
to those obtained with long characteristics codes: although up- 
wind intensities do not need interpolation along the beam, dif- 
fusion affects the local flux divergences when transfered from 
the slanted long characteristics grid back to the hydrodynamical 
grid. 

The basic mesh on which radiative transfer is computed is 
imposed by the MHD solver This is usually not critical in the 
optically thin upper atmosphere and the optically thick inte- 
rior, where radiative transfer is simple and may even be over- 
resolved. The opposite is the case in the transition region around 
the optical surface, where opacities drop rapidly due to their 
stron g temper ature dependence and cause a runaway cooling ef- 
fect JStein & N ordlund 199 q). For a s olar simulation, ID tests 
performed by Nordlund & SteinI (11991 ) indicate that a vertical 
spacing of < 10 km is desirable at this atmospheric height. Using 
a non-linear vertical grid with the finest resolution around the 
surface, this is easily achievable in 3D for modern MPI-based 
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Fig. 3. Numerical difFusion of a searchlight beam with rectangular cross-section using linear (upper right), local cubic (lower left) 
and local cubic monotonic (lower right) interpolation, compared to the exact solution (upper left). 



domain decomposed radiative hydrodynamics codes. However, 
for large coronal simulations or in the case of giant stars, where 
the spatial scales needed to resolve hydrodynamics and radiation 
transport exhibit much larger disparity than in the Sun, finding 
the optimal grid leads to a conflict. Besides the larger simula- 
tion size, too small length intervals Ax drastically increase the 
stiffness of the hydrodynamical equations, where the stability- 
limited time steps of the transport and diffusion terms scale 
with Ajc and Ajc^, respectively, and quickly become exceedingly 
small. In extreme cases, both effects may increase computation 
times of a model beyond tractabiUty. 

A fully adaptive mesh for computing radiative transfer would 
yield optimal results without affecting the stiffness of the equa- 
tions, but is difficult to realize in a characteristics method. We 
achieve partial adaptivity by inserting horizontal layers in the 
hydrodynamical mesh for the radiative transfer computation, re- 
ducing optical path lengths without reducing the time steps. The 
refinement is based on the maximum vertical gradient of the 
Rosseland mean opacity in each layer and reassessed in regular 
intervals. While inserting additional layers slows down conver- 
gence of the Gauss-Seidel method (see Sect. I3.2l l. this is again 
overcome by storing the source function from the previous time 
step. 



3.4. Numerical flux divergences 

Having established a method for numerically computing radia- 
tive transfer with coherent scattering in a decomposed simula- 
tion domain, we now need to obtain flux divergences V • f , a 
derivative of the radiation field. 

The right-hand side of Eq. (|7]i involves only local quantities 
that are defined on the cell centers of the hydrodynamical mesh, 
where Qi-^^ is eventually needed, and therefore seems a natural 
choice. The expression ;i^(7-5) is numerically stable in the opti- 
cally thin regime, where round-off errors of a possibly vanishing 
difference between J and S are attenuated by the exponential 
outward decrease of the opacity ;^'. At the same time, x amplifies 
round-off errors of (J - S) beneath the optical surface, where 
the radiation field thermalizes (J ^ B, also in the scattering case 
since e > 0): the flux divergence again vanishes, but the finite 
machine precision prevents complete cancellation of the terms. 

It is possible to stabilize a short characteristics solver in the 
whole simulation domain by subtracting So from the discretized 
formal solution (Eq. (fT4l i). which yields the modified integra- 
tion constant ^'^ - To - 1- Using this equation, one obtains 
a monochromatic Qmd.Aix, h) along each ray. We note, how- 
ever, that this leads to a deviation between the radiative energy 
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Jj Qmdjix, n)dD.dV that is emitted by the gas in the simulation 
volume V per time unit, and the outgoing radiative flux com- 
puted from the specific intensities at the surface: the expressions 
are not equivalent anymore in their discretized form, and numer- 
ical errors affect the two values in a different way. 

The discretized flux divergence V • F on the left-hand side 
of Eq. (|7| using finite difference quotients is stable in the opti- 
cally thick regime, but its accuracy deteriorates outward: round- 
off errors quickly become significant, as the internal energy per 
gas volume decr eases exponentially (see also the discussion in 
iBruls et alJ[T999l) . 

Adopting the a pproach presented in iBruls et all d 19991) and 
IVogler et alJ (|2005l) . we combine both expressions through ex- 
ponential bridging in each vertical column of the simulation do- 
main as a function of bin optical depth to benefit from their re- 
spective advantages. We slightly reduce the transition range be- 
tween the regimes by a squared exponent, resulting in the ex- 
pression: 



a-ad = .-^^^-^^e4d + (i-^-^^^"^G.' 



F 

■ad' 



(22) 



where tq = 0.1, Q^|.^^ - AnxiJ-S) and Q^^^ = -V • F, represent- 
ing the two sides of Eq. ^. The total radiative energy computed 
with this expression delivers a consistent surface flux, since Q^^^ 
contributes m ost of the radiative h eating. 

Following I Vogler et alJ (l2005l) . we compute radiative trans- 
fer on cell corners to improve the accuracy of g^^. Radiative 
fluxes F are averaged over cell corners surrounding each face be- 
fore computing difference quotients, while Qjf^^ is averaged over 
all eight cell corners surrounding each grid point. Both expres- 
sions use exactly the same stencil and exhibit very good agree- 
ment around the threshold optical depth tq in our solar-type sim- 
ulation. 

Flux divergences are computed only on the hydrodynamical 
grid. Additional layers that are possibly inserted by the radiative 
transfer solver just serve to stabilize the computation and may 
simply be omitted when computing Q^d, since conservation of 
the radiative energy flux through the hydrodynamical cell sur- 
faces must hold. 



4. Absorption and scattering opacity sources in the 
Sun 

A complete description of radiative transfer in stellar atmo- 
spheres requires a detailed wavelength-resolved treatment of 
numerous radiative absorption and emission processes, colli- 
sions with neutral atoms, electrons and ions in the plasma, as 
well as an evaluation of the feedback of the radiation field 
on the level populations of the interacting particles. The com- 
plexity of the resulting problem vastly exceeds current com- 
putational resources. We therefore restrict all of the underly- 
ing thermodynamical plasma states to LTE, neglecting the ef- 
fects of radiation on the excitation and ionization of atoms and 
photo-dissociation of molecules. The cross-sections and level 
populations needed for the absorption and scattering coeffi- 
cients then depend only on the gas density p and the temper- 
ature T. Microscopic plasma thermodynamics is treated with 
the Mihalas-Hummer-Dappen equation of state (EOS) for stel- 
lar envelopes (Hummer & Mihalas 1988: iMihalas etalj [19881 
iDappen et anil988t iMihalas et alJll99 0) and used in tabulated 
form. The solar chemical composition for the 15 elements in- 
cluded in the EOS and for the opacities is taken from the abun- 
dances of lAsplund et alJ (l2005h . 



4.1. Continuum opacity 

The most important continuous opacity sources are various tran- 
sitions of hydrogen atoms, their ions and molecules. The H" ion- 
ization opacity dominates the solar continuum around the optical 
surface in the visual band; the large temperature sensitivity of 
the weakly bound second electron in the hydrogen atom causes 
runaway radiation cooling and the strong tem perature gradient 
found at the top of the granules in the Sun (S tein & Nordlundl 
199 8). Most solar continuum photons originate from this very 
thin layer Among many other processes, photoionization of met- 
als contributes significantly to the continuous opacity at shorter 
wavelengths. Table lD.fl gives an overview of all sources consid- 
ered in this work; our data is mostly i dentical to those used i n 
the latest MARCS models (see Table 1 inlGustafsson et alj|2008b, 
but includes additional bound-free data from the Opacity Project 
and the Iron Project (see Trampedach et al., in prep., for further 
details). We also include opacities of the second ionization stage 
for many metals, allowing 3D models to extend deeper into the 
convection zone than their ID counterparts, which is a require- 
ment for correctly simulating surface granulation. 

The upper panel in Fig. |4] shows the wavelength and depth 
dependence of the continuum photon destruction probabilities 
e^ for the mean stratification of our 3D model, including all 
continuous absorption and scattering opacity sources consid- 
ered here. Continuum scattering has a significant contribution 
mostly above the surface, photons thermalize beneath at al- 
most all wavelengths. Note that the narrow features at the short- 
wavelength end are the scattering resonances of the Lyman se- 
ries; Lyman lines are nevertheless treated as true absorbers if 
line scattering is not included in the simulations. The Rayleigh 
scattering tail of HI contributes mostly to the UV continuum 
opacity in the upper solar photosphere due to its comparatively 
small cross-section and strong wavelength dependence (cta ~ 
/I"'*). The importance of elastic scattering on neutral hydrogen 
is outweighed by thermalizing processes closer to the surface 
and at short wavelengths. Electron scattering is wavelength- 
independent in the spectral range considered here, and becomes 
significant in the upper photosphere, where metals are the most 
important electron donors. It is mostly notable red-ward of the 
1 .644 yum edge of H" bound-free, before H" free-free absorption 
takes over Rayleigh scattering on He I atoms only gives minor 
contributions to the UV continuum opacity in the upper pho- 
tosphere. The scattering opacity of H2 molecules is negligible. 
Rayleigh and electron scattering are treated as isotropic, neglect- 
ing their weak (1 -H cos^ 6) anisotropy, where 6 is the scatterin g 
angle away from the incident direction (see, e.g.. lMihalasl 19781) . 

Between 5000 A and 1 .644 yum, the strong H" bound-free ab- 
sorption opacity thermalizes the photons. Its dominance slightly 
decreases in the cool outermost layers owing to the lack of free 
electrons to form the ion. 



4.2. Line opacity 

Spectral line absorption and scattering are important processes 
which dictate the near-radiative equilibrium found in the solar 
photosphere. The heating/cooling effect of this line-blanketing 
forces the flatness of the observed temperature gradient, bal- 
ancing the adiabatic dynamical gradient; see the discussion 
in Sect. 15.41 Spectral lines are particularly significant opacity 
sources at short wavelengths where many radiative bound-bound 
transitions of metals lie. 

We obtain line opacities from extensive opacity sampling ta- 
bles provided by B. Plez (2008, priv. comm.) as part of the MARCS 
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Fig. 4. Wavelength and depth dependence of the continuum photon destruction probabilities e^ (upper panel), the Ivan Regemortej 
(Il962h line photon destruction probabilities e', (center panel) and the total photon destruction probabilities q (lower panel) for the 
mean solar-type stratification. The zero point on the depth axis marks the stellar surface. 



collabora tion. The data are based o n VALD with some modifica- 
tions; see iGustafsson et all ( l2008l) for further details. The origi- 
nal line data combine scattering and absorption contributions in 
a total opacity, which is sampled with ~ 100000 wavelengths 
and tabulated for a range of temperatures and pressures. The ta- 
bles assume Saha ionization equilibrium and Boltzmann level 
populations to obtain the absorber density fractions. Departures 
from LTE, e.g. through radiative ionization, are neglected. 

Following ISkartlierJ ( 120001) . we estimate the importance of 
scattering in line transitions by computing a photon destruc- 
tion probability e', for every line opacity sample, using the 
[van Regemorter (1962) formula (see Appendix ICl). We assume 
all scattering atoms to be neutral, accoun ting for the large con- 
tribution of Fel to the line-blanketing ( Anderson! 1 19891) , and 
all tra nsitions to be permitted , in which case the assumptions 
of the Ivan Re gemorteij (Il962h formula yield reasonable esti- 
mates. Only electrons are taken into account for collisional de- 
excitation. The estimated photon destruction probability e^ is 
then a function of wavelength, temperature and electron pres- 
sure, and independent of the actual transition. It may therefore 
also be applied in cases where the line opacity sample includes 
several transitions (see the discussion in AppendixICli. Line tran- 
sitions are treated as independent two-level processes without 
taking the coupling of the respective level populations into ac- 
count, which is a reasonable assumption for resonance lines. 

The center and lower panels in Fig. |4] show the wavelength 
and depth dependence of the estimated e'^ of spectral lines and 
the total photon destruction probabilities e^, including all con- 
sidered continuous and line processes. It is clear that collisional 
de-excitation dominates beneath the surface and at the longest 
wavelengths. Resonant line scattering becomes important to- 
wards optical and shorter wavelengths at increasing depth. 

With the exception of very strong lines, line scattering is gen- 
erally not coherent due to the Doppler shifts in the moving gas, 
which are not accounted for in our calculations. The two-level 



approximation probably gives a reasonably realistic picture of 
strong permitted lines, but departures from the LTE populations 
of the atomic levels are still neglected. The important Fe I opac- 
ity deviatesjromjhe LTE estimate in higher layers (see Fig. 7 
in lShort & Hauschi ldt 2005), thereby affecting the overall mag- 
nitude of the line-blanketing in these regions. Moreover, the ac- 
curacy of the opacity sampling method itself deteriorates out- 
wards, where fewer and fewer lines contribute to the opacity. 
The van Regemorter approximation assumes resonant line scat- 
tering and consequently produces poorer estimates for all non- 
resonant lines. In summary, we should expect to obtain an order- 
of-magnitude estimate for the effects of scattering on the atmo- 
spheric structure. A more detailed picture requires a full treat- 
ment of the departures from LTE level populations and velocity 
fields, which is still out of reach for time-dependent 3D simula- 
tions. 



5. The effects of scattering on the photospheric 
temperature structure of a solar-type star 

5. 1 . The 3D hydrodynamical surface convection model 

To investigate the effects of scattering on the atmosphere of a 
solar-type star, we conduct time-dependent radiative hydrody- 
namical simulations of the quiet surface, neglecting the effects of 
magnetic fields. We solve the fully compressible Navier-Stokes 
equations, the mass conservation equation and the energy equa- 
tion, along with the time-independent radiative transfer equation 
(Eq. dUi); see, e.g.. iStein & Nordlundl(ll998l) and lNordlund et aP 
(.2009) for further details. Our 240 x 240 x 226 model covers 
a horizontal area of 6 Mm x 6 Mm at a constant resolution of 
25 km, and extends approximately 700km above and 2.8 Mm 
below the surface. The vertical resolution reaches 7 km around 
the radiative cooling peak and decreases in the optically thick 
and thin parts of the simulation; radiative transfer is thus re- 
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Fig. 5. Horizontal average heating rate per unit mass around the 
stellar surface at an arbitrary time step of the simulation. Boxes 
show (^lad) when computed on the hydrodynamical grid; the 
vertical resolution reaches 7 km around the peak. The solid line 
shows the result after inserting two additional horizontal layers 
in each hydrodynamical cell. The upper panel gives the average 
deviation between the two cases in the same units. 



solved well enough that only ~ 3 % of the rays would be affected 
by overshoots (see Sect. 13. lb . We test the accuracy of the verti- 
cal resolution using the adaptive refinement, inserting two extra 
layers before each computation of radiative transfer Local dif- 
ferences between the two calculations reach ~ 3T0"'ergg"' s"', 
owing to the strong sensitivity of the heating rate per unit mass, 
^rad = Qmilp, to the local temperature gradients in the highly 
inhomogeneous granulation flow. On the average, however, the 
change in radiative flux divergence is negligible (see the upper 
panel of Fig. |5]), and the radiation field is well resolved on the 
hydrodynamical grid. Note the difference between the magni- 
tude of the cooling peaks in Fig.|5]and Fig.|6| the ID calculation 
is based on the mean structure; in the 3D case, the average over 
each depth layer in the 3D box is taken and thus includes lateral 
inhomogeneities produced by the granulation flow. 

Horizontal boundaries are periodic to mimic an infinitely ex- 
tended atmosphere, vertical boundaries at the top and bottom of 
the simulation box are open to minimize the interference with 
the granulation flow. Mass conservation is ensured at the bottom 
by keeping the gas pressure constant; the underlying convection 
zone is mimicked by setting the entropy of the inflowing gas. 
The upper atmosphere is stabilized by setting internal energies 
to a slowly evolving average at the top. 

We approximate the wavelength integral (Eq. ([8])) with 
12 opacity bins to account for the depth-dependence and 
wavelength-dependence of the absorption and scattering coeffi- 
cients. The simulation box extends far into the optically thin at- 
mosphere with (tsooo) ~ 10"^, where irradiation I^ from above 
is negligible. Rosseland optical depths at the bottom typically 
reach (tross) ~ 10^, where radiative transfer is entirely diffusive 
and the radiation field is completely thermalized. We therefore 



set the diffusion approximation I^^^ - Bbot + dB/dr for all ingo- 
ing intensities at the bottom. 

The three simulations discussed in Sect. 15.31 have mean ef- 
fective temperatures Tgff between 5804 K and 5811 K with av- 
erage temporal fluctuations of about 1 3 K; they are thus slightly 
hotter than the Sun. For our purposes, there is no need to ex- 
actly reproduce the solar Teff . The simulations yield time-series 
of snapshots spanning ~ 1 h of stellar time each, covering sev- 
eral granule lifetimes (t ~ lOmin) and several periods of the 
dominant p-mode (f ~ 5 min). Our simulation box covers about 
10 granules with typical sizes of the order of ~ 1 Mm, allow- 
ing us to obtain a statistically meaningful sample of the surface 
flow in terms of the ergodic hypothesis. The model without scat- 
tering was computed with a coarser radiation time step of 0.2 s, 
keeping the radiation field constant during the intermediate hy- 
drodynamical calculations. The slow evolution of the flow field 
and the locality of the Planck source function allow such reduc- 
tion of the computation time in very good approximation. 

5.2. Scattering in tiie 1D mean stratification 

We first test the importance of scattering in the ID mean strat- 
ification of our 3D model (the S = B case, see Sect. 15.31 ) by 
comparing the wavelength-integrated ^lad, using the full opacity- 
sampled spectrum. Radiative transfer was computed in ID using 
a direct block matrix Feautrier-type solver with coherent scat- 
tering (for a detailed description see, e.g., Rutten 2003) and 4th 
order Radau quadrature for the integral over the polar angle. The 
left-hand side of Fig.|6]shows ^rad without scattering and S - B, 
with continuum scattering only, and with both continuum and 
line scattering (lower panel), as well as the deviations from the 
first case (upper panel). 

Continuum scattering seems to have very little impact on 
Imd for the given mean structure; the cooling is slightly stronger 
near the surface. This behavior is expected from the mostly large 
photon destruction probabilities e^ shown in the upper panel of 

Fig. a 

The differences are slightly larger when scattering is in- 
cluded in the line-blanketing: the small heating bump, where 
cool uprising ga s is heated from beneath by hot granules (see 
the discussion in IStein & Nordlundll998h . and the cooling peak 
beneath the surface both slightly weaken, since the fraction of 
scattered photons in the line-blanketing does not contribute to 
heat exchange (cf. the right-hand side of Eq. Q). The upper at- 
mosphere, however, now shows slight heating of the mean struc- 
ture. 

We repeat the same test with the binned opacities, comput- 
ing ID radiative transfer with and without scattering for 12 mean 
opacities, photon destruction probabilities and bin-integrated 
Planck functions. The right panels of Fig. |6] compare again the 
three different cases. The binning has been optimized for match- 
ing sampled and binned ^lad in the S = B case (solid lines in 
the lower panels of Fig. |6]l. The continuum scattering calcula- 
tion with opacity bins underestimates the cooling beneath the 
surface. The disparity increases further when line scattering is 
included; the relative deviations reach 7.5 % in the cooling peak 
(dot-dashed lines in Fig.|6]l. However, the overall impact of scat- 
tering radiative transfer on the temperature structure of the 3D at- 
mosphere above T5000 ^ 10"^ is small (see Sect. l5.3l and Fig.|7]l, 
the same binning setup was therefore adopted for all three sim- 
ulations. Higher up in the atmosphere, at tsoqo < lO""', opacity 
binned radiative transfer shows slightly stronger heating of the 
gas. 



12 



W. Hayek et al.: Radiative transfer witii scattering for domain-decomposed 3D MHD simulations of cool stellar atmospheres 



en 



0.8 
0.6 
0.4 
0.2 
0.0 
-0.2 

0.0 
-2.0 



2 -4.0 



-6.0 



-8.0 



10" 



- 




- 


L 


/ 'x 


] 


L 




■ 






: 












. 




— ^ 


r-: 




^ 




— No scot. / 


■ 


- 


- - Cont. scot. / 


- 




-■- Cont. + Line scot. / 


■ 




V 

1 1 1 1 1 1 


■ 



0.8 
0.6 

I 0.4 

< 0.2 

0.0 

-0.2 

0.0 

en 

'c" -2.0 
a^ 

o 

£ -4.0 

.E 

T> 
O 

^ -6.0 



10"" 10" 



10"^ 



10" 
sooo 



10° 10' 



10=^ 



-8.0 



10" 






— No scot. 

- - Cont. scot. 
-■- Cont. + Line scot. 




10"" 10" 



10" 



10"' 10° 



10' 



10=^ 



Fig. 6. Left: Heating rates ^lad per unit mass as a function of monochromatic optical depth at 5000 A, computed on the ID mean 
structure with full opacity sampling for three cases: without scattering (solid line), with continuum scattering (dashed line), and with 
continuum and line scattering (dot-dashed line). The upper panel shows the deviations of the latter two cases from the computation 
without scattering (same axis units as in the lower panel). Right: Same computation, but using mean opacities and scattering albedos 
in 12 bins for the radiative transfer computations. 



5.3. Scattering in the mean 3D model 

In order to assess the effects of continuum and line scattering, we 
perform three independent simulation runs: the first one treats 
radiation without scattering by adding all scattering opacity to 
the absorption opacity and assuming a Planck source function 
S = B. The second one includes continuum scattering in the 
source function and only adds line scattering opacity to the ab- 
sorption opacity, and the third one includes scattering both in 
the continuum and in the line-blanketing. All three time series 
start from the same initial snapshot and span the exact same 
amount of simulation time. Snapshots are taken at regular in- 
tervals of Afsini = 10 s. We consider time steps at fsim > 8min 
after the initial snapshot to allow the atmosphere to adjust to any 
changes in the radiative heating rates. Exploiting the tight cor- 
re lation between gas temp erature T and vertical optical depth 
T dStein & Nordlundll 19981) . we interpolate the 3D temperature 
cube at each time step of the series onto surfaces with the same 
optical depth, using a reference r-scale at 5000 A. We then com- 
pute the average temperature of each surface in the 3D cube, 
which yields a ID mean temperature profile for every snapshot. 
These profiles are finally averaged over time, and we obtain a 
very robust characteristic T-t relation. 

Figure [T] compares the resulting horizontal and temporal 
mean temperature profiles. The simulations without scattering 
and with continuum scattering have practically identical strat- 
ifications, as expected from the continuum photon destruction 
probabilities e^ (Fig.|4| and the ID test presented in the previ- 
ous section; continuum scattering is therefore insignificant for 
the atmospheric stratification in solar-type stars. 

The effects of scattering on line-blanketing in and below the 
photosphere are also rather weak (dot-dashed line in Fig.|7j- The 
gas temperatures above T5000 Z 10"^ deviate up to 40 K from the 



stratification without scattering, resulting in a slightly steeper 
temperature gradient around the surface (tsqoo = 1). Since our 
adopted binning setup overestimates the deviations for the ID 
mean structure (right-hand side of Fig. |6]), the impact of line 
scattering is probably even smaller at tsooo <; 10"^. The temper- 
ature structure in the lower photosphere is thus hardly affected 
by scattering. The opposite is the case in the high photosphere 
and above (T5000 ^ 10""*), where we observe temperatures that 
are about 350 K lower, resulting in a significantly steeper mean 
gradient. 



5.4. Comparison of the Wand 3D calculations and with 
other model atmospheres 

The effects of line scattering on the temperature structure of the 
3D model seem to be opposite of ID hydrostatic models in radia- 
tive equilibrium, where heating of the highest layers rather than 
cooling is observed. Indeed, the ID calculations on the mean 3D 
atmosphere exhibit slight heating in this region when scattering 
is included (Fig. |6]l. The temperature gradient would therefore 
become shallower if the ID calculations were iterated under the 
assumption of radiative equilibrium (see, e.g., the discussion in 
iRuttenll2003l) . 

The total radiative flux divergence includes several com- 
ponents: hot radiation from deeper layers at short wavelengths 
dominates the heating of the gas; the steep outward dB,\ldT 
gradient causes a positive growing (J - S) split. The effect de- 
clines in higher layers due to the rapidly decreasing opacity (cf. 
Eq. dTjl). Strong LTE lines may heat or cool the higher atmo- 
sphere (since J x B in deeper parts), depending on the spec- 
tral region and local temperature gradient, which determine the 
sign of the (J - S) split. Including coherent scattering in line- 
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blanketing effectively reduces both radiative heating and cooling 
in high layers through the outwards decreasing e' (see Fig.[T]). As 
a consequence, strong resonance lines become unimportant for 
the temperature structure in high layers, and radiative heating at 
shorter wavelengths decreases. 

In the ID mean atmosphere, scattering-weakened line cool- 
ing shifts the total q^a^ slightly towards positive values. The be- 
havior of the 3D case can be understood by considering the 
dynamical nature of our 3 D models. Following a derivation in 
iMihalas & Mihalail (Il984l) . we insert the continuity equation 



Dp 

-^+pV -u =0, 
Dt ^ 



(23) 



where p is the gas density, D/Dt is the material derivative and u 
is the gas velocity, into the energy equation. 



De P 
Dt p 



(24) 



where e is the internal energy per unit mass, P is the gas pressure, 
and qaA is the radiative heating rate per unit mass; we omit the 
viscous dissipation term for simplicity. The resulting expression. 



De_ 
'Dt 



P Dp 



-=- 7-r— - ?iad, 



(25) 



is the first law of thermodynamics. An upflowing (downflowing) 
gas parcel cools (heats) through expansion (compression) repre- 
sented by the Dp/Dt term in Eq. (IZST i, and is exposed to radiative 
heating through the qi-^d term. Equation dZSl l is equivalent to the 
expression 



Ds 

T — 

Dt 



qtai, 



(26) 



where T is the gas temperature and s is the entropy per unit mass, 
and it is immediately clear that gas motion is adiabatic when 
^rad — * 0. In the photosphere of the 3D simulation, temperatures 
are not affected by scattering. In the upper atmosphere, below 
T5000 ~ 10"'*, scattering strongly reduces the line-blanketing. 
Small or vanishing heating rates ^lad cause the temperature strat- 
ification to steepen towards an adiabatic gradient. 

Figure |8] compares the time evolution of radiative heating 
rates at three optical depths in the atmosphere, averaged over sur- 
faces of constant optical depth to approximately account for ver- 
tical gas motion. The plot shows a sequence of snapshots taken 
at regular simulation time intervals of 10 s; the zero point on the 
abscissa is arbitrary. At T5000 = lO""* and T5000 - 10"'' (lower 
and center panels), the continuum scattering case (circles) and 
the continuum and line scattering case (diamonds) exhibit simi- 
lar positive heating rates on the average (dashed and dot-dashed 
lines) and thus similar average temperatures (Fig.|7]). Line scat- 
tering radiative transfer produces slightly stronger mean heat- 
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ing at T5000 - lO"-', but fluctuates with lower amplitude. At 
T5000 = 10"^, ^rad practically vanishes on the time average in the 
line scattering case, but there is still significant radiative heat- 
ing with line scattering as true absorption. Note the dynamical 
variation of the sequences: contrary to ID hydrostatic models, 
where the radiation field is time-independent by definition, the 
evolution of the 3D simulations produces fluctuating radiative 
heating. 



IWedemever et al.l (l2004t) presented 3D radiation- 
hydrodynamical simulations of the solar atmosphere that 
include a chromosphere, using radiative transfer without scat- 
tering and solving the equation only for the Rosseland mean 
opacity to suppress radiative cooling by strong LTE lines. They 
found an increasing asymmetry of the gas temperature distribu- 
tion with increasing height above the surface, and a bifurcation 
in the chromosphere. IWedemever et al.l (l2004l) further observed 
that treating strong spectral lines as true absorption with the 
opacity binning method reduces the amplitude of temperature 
fluctuations, which are caused by outward propagating acoustic 
waves, resulti ng in unrealistica lly low maximum temperatures 
in high layers. ISkartlieiJ (l2000h investigated scattering radiative 
transfer in the chromosphere, comparing radiative heating 
with and without scattering, and came to the conclusion that 
including line scattering reduces this damping effect of LTE 
lines. 

Our simulations do not include a chromosphere; the inter- 
nal energy at the top boundary is set to a slowly evolving mean 
instead. In the line scattering case, where radiative transfer has 
only weak influence on the gas, the temperature gradient is sen- 
sitive to this boundary condition and thus not well-constrained. 
However, this does not compromise our conclusions, since the 
boundary is free to adapt to any upward or downward shift in the 
mean energies of the gas beneath. 

Figure |9] shows temperature distributions of the simulations 
with continuum scattering and with continuum and line scatter- 
ing at three different heights above the surface. Our si mulations 
do not reac h the same geome t rical h eights as those of lSkartlienl 
(200(t) and ' Wedemever et al.l (l2004l) . and we use a more real- 
istic radiative transfer treatment with 12 opacity bins. We find 
a similarly growing asymmetry in the temperature distribution 
of the hne scatterin g simulation in the outer layers (cf. Fig. 7 in 
IWedemever et al.l20 04). Treating strong lines as absorbers shifts 
the mean temperature upward and removes the high temperature 
tai l of the di s tribut ion, i n qualitative agreement w ith the findings 
of lSkartlienl (|2000|) and lWedemever et all ( |2004 . 

Figure ITOlshows horizontal and temporal averages of the rel- 
ative temperature fluctuations, which we define as 



Ar™, ^J{{T - <r))2) 



(T) 
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(27) 



W. Hayek et al.: Radiative transfer witii scattering for domain-decomposed 3D MHD simulations of cool stellar atmospheres 



15 



in every geometrical depth layer (cf. Eq. 2 and Fig. 9 in 
IWedemever et alJ2004l) . The comparison between the cases with 
continuum scattering and with continuum and line scattering 
confirms the damping of temperature fluctuations through line 
absorption. Note the decreasing ATnns at the top of the simula- 
tion, which is induced by the hydrodynamical boundary condi- 
tions. 
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We conclude that line scattering is an important ingredient 
for model atmospheres of solar-type stars that include a chromo- 
sphere; while gray radiative transfer reduces damping through 
strong LTE lines, it cannot produce a reaUstic temperature struc- 

turC; 

Anderson (|l989h presented simplified, and 

IShort & Hauschildll (l2005l) presented full ID non-LTE 
line-blanketing calculations, respectively, for hydrostatic 
model atmospheres of solar- type stars. The departures of 



line-blanketing from LTE through iron-group elements heat 
up the atmosphere in the height range 10"^ < T5000 ^ 10"^. 
Our 3D model predicts a predominant temperature decrease as 
we discussed above. However, it is not clear how departures 
from LTE in the absorber populations through the ionization 
balance etc. would affect the atmospheric structure in our 3D 
simulations, making a direct comparison with the ID models 
difficult. 

Doppler shifts may have a significant influence on line ab- 
sorption in higher l ayers, which e xpose line cores to hot radia- 
tion from deeper in. Vogler et al. (2004) estimated the effects to 
be insignificant in the photosphere, but their work was based on 
ID tests. The large scattering albedo of strong resonance lines, 
however, should reduce the impact of Doppler shifts higher up. 

6. Conclusions 

We presented a 3D radiative transfer method with coherent scat- 
tering for time-dependent (M)HD simulations of stellar atmo- 
spheres with the new BIFROST code (Gudiksen et al., in prep.). 
The simulations are parallelized through domain-decomposition 
to take advantage of large-scale computer clusters. The solver is 
based on short characteristics and the Gauss-Seidel scheme for 
an iterative computation of the radiation field and the radiative 
flux divergence in the whole simulation domain. We use mono- 



tonic interpolation to reduce the numerical diffusion effect of 
short characteristics and represent the source function integral 
with Bezier polynomials to suppress interpolation overshoots. A 
partial grid refinement scheme is included to improve the reso- 
lution of the radiative transfer computation where strong verti- 
cal opacity gradients occur The wavelength integral is treated in 
the opacity binning approximation, using 12 bins that divide the 
opacity spectrum by formation height and wavelength. 

The effects of coherent scattering on the temperature struc- 
ture of a solar-type star are investigated with 3D time-dependent 
hydrodynamical simulations of magnetically quiet surface con- 
vection, including Rayleigh scattering and electron scattering 
in the continuum and estimated line scattering using the van 
Regemorter formula. While continuum scattering processes are 
not important for the mean temperature stratification, we find 
lower temperatures in the upper atmosphere when scattering is 
included in the line-blanketing. 3D radiative-hydrodynamical at- 
mospheres thus show the opposite behavior of ID hydrostatic 
atmospheres in radiative equilibrium, where scattering in strong 
lines effectively heats the outer layers. 

3D LTE models of solar surface convection have been very 
successful at reproducing various observational tests, and our 
results indicate that the solar photosphere is indeed well rep- 
resented when scattering is not included in radiative transfer It 
therefore seems that a refined treatment of the line-blanketing 
through, e.g., opacity distribution functions or opacity sampling 
will be the next significant step to improve the realism of 3D 
radiative-hydrodynamical model atmospheres. Scattering radia- 
tive transfer is nevertheless an important ingredient of consistent 
3D MHD models of the solar chromosphere, transition region 
and corona. 

While it is not unexpected to see only small differences in 
the photospheres of solar-type stars when scattering is taken into 
account, this is likely to change for the much less dense atmo- 
spheres of giants, where the importance of Rayleigh scattering 
increases. The case of metal-poor giants is particularly interest- 
ing in that respect, owing to their significance for understanding 
galactic chemical evolution and the origin of the elements. 
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Appendix A: Bezier interpolation of source 
functions and opacities 

The discrete foimal solution (Eq. )14t ) of the radiative transfer equation (Eq. J4)) 
requires interpolating the source function S(t) along the short characteristic. 
While linear interpolation never overshoots, its accuracy is not sufficient in opti- 
cally thick media, since the discrete expression is not equivalent to the diffusion 
approximation. Second-order interpolation significantly improves the accuracy, 
but suifers from strong overshoots where AS /At gradients change rapidly be- 
tween the upwind and downwind halves of the characteristic. In extreme cases, 
this can even destabilize the solver and produce spikes in the local flux diver- 
gences. 

Bezier-type interpolation techniques allow for a direct detection and sup- 
pression of such over shoots by virt ue of a control point in the polynomial which 
shapes its curve (see iAuenl2003l) . A second-order Bezier polynomial may be 
written in the parameterized form 



S(t) = Su{l-tf +S,/ + 2SM ■ 



0, 



(A.l) 



where Su and So are the source functions at the upwind end and the center point 
of the characteristic, between which interpolation is needed, t = (r-T„)/(To -Tu) 
is the curve parameter, and Sc is the control point. A Bezier curve is always 
bounded by the convex hull of the three points Su, Sc and 5o. Using the abbre- 
viations 

An, =To-ri,; ATd=rd-To (A.2) 

for the optical depths along the characteristic and choosing the control point 

At„ 



Sc 



Lis' 







At„ 



Atj So -Su 



At,, 



-So 



2 \ Atu + ATd Atu Atu + Atj Atj 



yields second-order interpolation of S , which now also depends on the source 
function Si at the downwind end. Introducing the three functions 

Uo(t) = l-e-' 
Um = t-Uo(t) 
U2(t) = t^-2Ui{t), 

and evaluating the integral of the Bezier polynomial results in the familiar 
second-order integration constants for Eq. J14t , 



{;2(AT„)-(ATd-^2ATu)t/l(ATu) 
At„(Atu +■ Atj) 
(At„ + Atj){/i(At„)-[/2(AtJ 

AT„ATd 

{/2(At„)-At„J/|(At„) 



>F„ = Uo^^r^) + 
'Po = 
I'd = 



ATd(ATu + ATd) 



(cf. Eq. 8 and 9 in lKunasz & Aueall988h . If the source functions Su, 5o and 
Sd have an extremum at So, choosing S^ = Sq enforces Sq = 0, yielding the 
constants 



T, = f/o(AT„) + 



(72(At„)-2At„[/i(At„) 



1-0 



At; 
2At„J/i(At„)-C/2(At„) 



At;^ 



I'd = 0. 
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Overshoots are avoided by limiting 5 c to the range of the data points: 
min(Su,So) < Sc < max(Su,So)- If Sc lies outside these limits, choosing 
Sc = Sa results in the constants 



"Vu 



"Vo 



'Frf 



f/o(Ar„) - 
[/2(Ar„) 



t/2(Aru) 
Arl 



Arr. 



0. 



Note that, contraiy to the first two cases, suppressing such overshoots leads to 
discontinuous left-hand and right-hand derivatives at 5 o ■ 

Optical depths are computed in an analogue fashion to avoid negative re- 
sults. A second-order Bezier polynomial ^(.v) interpolates opacities over the path 
length As along the ray; integration over .s yields the optical depth interval 



Ar 



r 

Jo 



X(cr)d<T -- 



As 



■Cfu+^O+^c), 



(A.3) 



where the control point Xc is selected according to the same criteria as discussed 
above for S c . 



Appendix B: Local cubic monotonic interpolation 

The radiative transfer solver uses local cubic interpolation for interpolating 
data from the hydrodynamical grid onto the characteristics grid. This choice of 
method improves the accuracy compared to linear interpolation, while ensuring 
local control of the interpolating polynomial to reduce artifacts. In addition to 
being a one-pass algorithm, the method also exhibits good computational per- 
formance thi'ough its instruction-per-data ratio, which is well-suited for modern 
multi-core CPUs, where high computation speeds are contrasted with slow mem- 
ory access. 

2D interpolation is approximated by consecutive ID interpolation using a 
cubic polynomial 

f(t) = aP +hl- + ct + d, (B.l) 

with the curve parameter ? e [0,1]. The coefficients a, b, c and d depend on the 
adjacent data points f\ and /i and their derivatives f'^ and /,'. Inserting the data 
and reordering the terms, the polynomial may be rewritten in the form 

fit) = a(t)f\ + mil + y(t)f[ + 5(0/2 . (B-2) 

where the interpolation weights a, /S, y and S now depend on the pai'ameter (: 

a(t) = 2? - 3r^ + 1 



/?(/) = 3f^ - if" 

y{t) = it^ - 2r + t) Ax 



m 



■t^)Ax 



with the grid spacing Ax between the two data points. The shape of the curve 
is defined by the derivatives /[ and /^ . A natural choice is the mean of the left- 
handed and right-handed difference quotients f^ and f^ at both end points. An 
unweighted arithmetic mean leads to wiggles and overshoots where stro ng gra- 
dients appeal'. We therefore adopt a recipe bv lFritsch & ButlanJ 419841) . which 
uses a weighted harmonic mean 

/ f f 

' ikLB, f ' f > 

(l-a)f[+af^ JlJr 

I /l/r ^ 

with the weighting factor 



/' 



(B.3) 



1 + ■ 



AxR 



Aa'l + Axr 



(B.4) 



which depends on the grid spacing A.vl and A.itr on the left and right sides of 
the data point. The weighted harmonic mean biases /' towards the smaller of 
the two difference quotients f^^ and /^ where strong gradients occur, effectively 
suppressing overshoots. 

Quadratic interpolation uses only one of the two derivatives /[ and /^ , de- 
pending on the interpolation parameter /. The interpolation coefficients are 



t<^_: r>f : 

a{t) = 1-f ait) = (1 - 0^ 

m = t^ m = tn-t) 

rit) = til -t) Ax r(0 = 

Sit) = Sit) = tit - l)Ax 



(B.5) 



Appendix C: Line scattering with the van 
Regemorter formula 
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Fig. C.l. Dependence of the tabulated Gaunt factor integral 
P(AE;i , T) for co llisions of electrons with neutral atoms in the 
Ivan Regemorteij (11962.) formula on the transition energy AE^ 
and the gas temperature T. The dashed lines mark the boundaries 
for typical values of AE^IkT found in solar surface convection 
simulations (Sect. lSTt . 



The photon destruction probabilities in line transitions may be estimated 
using the semi-empirical van Regemorter 1 1962) formula to obtain electron col- 
lision rates, following the discussion in Skartlien (2000). Neglecting other con- 
tributions from, e.g., collisions with neutral hydrogen atoms, the de-excitation 
rate for electron collisions according to this formula is given by 



C2i~A^N^T-^I^A2iPiAE,,T), 



(C.l) 



where A is the transition wavelength, N^ is the electron density, T is the gas 
temperature, and /I21 is the Einstein coefficient for the corresponding sponta- 
neous radiative transition. The function PiAE^, T) abbreviates the velocity inte- 
gral over the empirically calibrated Gaunt factor of the scattered electron, and 
depends on the transition energy AEi and the gas temperature T. We adopt the 
tabulated data for neutral atoms of van Regemorter ( 1962), see Fig. lC.ll 
The photon destruction probability for a two-level atom is given by 



C21 



ka+o-a C21+A21+B21BA 



(C.2) 



where B21 is the rate for induced de-excitation, and B,i is the Planck function. 
Neglecting the induced de-excitation term, Eq. jC.2t simplifies to 



1 



I+A21/C21 



(C.3) 



e,i is independent of the actual transition after inserting the van Regemorter for- 
mula Eq. )C. IK and thus only a function of A, N^ and T. 

Line opacities in stellar spectra often combine contributions from many 
transitions at a given wavelength. The total monochromatic photon destruction 
probability of an opacity sample at wavelength A is given by the sum over all 
transitions. 



2^ 

^•^A.i 



(C.4) 



Inserting Eq. jC.3t , thus assuming the above mentioned approximations, yields 



i-'i^Aj 



(C.5) 



where the equality holds since f,i is independent of the actual transition i. The 
absorption and scattering contributions /?^ and o"'j to each opacity sample x\ are 
then isolated using ejj and added to the coefficients of the continuum processes 
(see TablelPHl. 
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Appendix D: Continuum opacity sources 



Table D.l. Continuum opacity sources 



Absorber and process Reference 



Hb-f 




Broad & Reinhardt (1976): 


;Wishart(1979) 


Hf-f 




Bell &Berrington (1987) 




HIb-f,f-f 




Karzas & Latter (1961) 




HI+HI 




Doyle (1968) 




Hl+Hel 




Gustafsson & Frommhold 


(2001) 


H2+HI 




Gustafsson & Frommhold 


(2003) 


Hz+Hel 


J0r£ensen et al. (2000) 




H2+H2 




Borvsowetal. (2001) 




H-f-f 




Bell (1980) 




H2 photo-dissociation 


Alhson Sl Dalaarno (19691 


1 


HJ b-f, f-f 




Standi (1994) 




He- f-f 




John (1995) 




He I b-f 




TOPbase' 




He I f-f 




Peach (1970) 




He II b-f 




TOPbase' 




c-f-f 




Bell etal. (1988^ 




CI b-f 




Nahar & Pradhan (1991) 




CI f-f 




Peach (1970) 




C II b-f 




Nahar (1995, 2002) 




C II f-f 




Peach (1970) 




cm b-f 




Nahar & Pradhan (1997) 




N-f-f 




Ramsbottom et al. (1992) 




NIb-f 




Nahar & Pradhan (1997) 




Nil b-f 




Nahar & Pradhan (1991) 




N III b-f 




Nahar & Pradhan (1997) 




o-f-f 




John (1975) 




01, on b-f 


Nahar (1998) 




III b-f 




Nahar & Pradhan (1994b) 




Nel, Nell, 


Ne III b-f 


TOPbase' 




Nal, Nail, 


Nam b-f 


TOPbase' 




MgI,MgI] 


[,MgIIIb-f 


TOPbase' 




AlI,AlII,AlIIIb-f 


TOPbase' 




Si I b-f 




Nahar & Pradhan (1993) 




Si II b-f 




Nahar (1995) 




Si III b-f 




TOPbase' 




SI b-f 




TOPbase' 




S II b-f 




Nahar (1995) 




S III b-f 




Nahar (2000) 




ArI,ArII, 


Ar III b-f 


TOPbase^ 




Cal, Call, 


Ca III b-f 


TOPbase' 




Felb-f 




Bautista(1997) 




Fell b-f 




Nahar & Pradhan (1994a) 




Fe III b-f 




Nahar (1996) 




Nillb-f 




Bautista(1999) 




CO- f-f 




.John (1975^ 




H2O- f-f 




.John (1975^ 




OH b-f 




Kurucz et al. (1987) 




CHb-f 




Kuruczetal. (1987) 





H I scattering 
H2 scattering 
e- scattering 
He I scattering 



Gavrila(1967) 

Victor & Dal gamo (1969) 

Thomson 

Langhoffetal. (1974) 



Contains Opacity Project data iSeaton et al.ll 1994h 



